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SECTION  A 


The  h-p  Version  of  the  Finite  Element  Method.  Theory,  Numerics 
and  Industrial  Experience 
I  Babuska,  University  of  Maryland,  U.S.A. 


The  talk  will  present  the  basic,  recent  results  of  the  h-p  version  of 
the  finite  element  method  and  the  applications.  Numerical  examples 
will  be  presented. 


The  Numerical  Solution  of  Initial-value  Problems  for  Integro-Pifferentiai 
Equations 

H  Brunner,  Memorial  University  of  Newfoundland,  Canada. 


In  this  talk  we  present  a  survey  of  recent  developments  in  the  numerical 
treatment  of  various  types  of  Volterra  integro-differential  equations  of 
order  r  i  1,  including  equations  where  the  r-th  derivative  of  the 
unknown  solution  occurs  also  in  the  kernel  of  the  integral  operator.  « 
The  discussion  will  focus  on  theoretical  and  numerical  aspects  of  spline 
collocation  methods  and  on  a  related  class  of  implicit  Runge-Kutta-Nystroni 
methods. 

The  talk  will  conclude  with  a  look  at  some  ongoing  work  regarding  the 
numerical  solution  of  partial  Volterra  integro-differential  equations 
of  parabolic  and  hyperbolic  type. 


The  Numerical  Solution  of  Nonlinear  Two-Point  Boundary  Value  Problems 
using  Iterated  Deferred  Corrections 
J  K  Cash,  Imperial  College,  London. 

Two  of  the  moat  popular  'global'  methods  for  solving  general  nonlinear 
two-point  boundary  value  problems  are  collocation  and  deferred  correction. 
For  first  order  systems  there  is  a  well  known  equivalence  between  certain 
collocation  and  Runge-Kutta  methods.  In  this  talk  we  examine  a  method 
which  combines  these  two  approaches  namely  Runge-Kutta  formulae  with 
deferred  corrections.  In  order  to  perform  the  deferred  correction 
stage  efficiently,  a  special  class  of  formulae,  known  as  Mono-Implicit 
Runge-Kutta  formulae,  is  derived.  The  results  of  some  extensive 
numerical  computations,  particularly  for  singular  perturbation  and 
turning  point  problems,  are  described. 
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Analysis  and  Computation  for  a  Model  for  Possible  Contamination  by 
Nuclear  Waste-Disposal  in  Porous  Media 
R  E  Ewing,  University  of  Wyoming,  U.S.A. 


The  transport  through  porous  media  of  contamination  generated  by 
nuclear  waste-disposal  cau  be  modelled  by  a  system  of  four  coupled 
nonlinear  partial  differential  equations.  The  total  fluid  flow 
is  described  via  an  elliptic  equation  heat  transport  by  a  parabolic 
equation,  and  the  miscible  displacement  of  the  dominant-  and  trace- 
species  by  transport-dominated  parabolic  equations.  Since  the  fluid 
pressure  appears  in  the  three  parabolic  equations  only  through  its 
velocity  field,  a  mixed  finite  element  method  is  used  to  approximate 
the  fluid  velocity  directly.  The  heat,  brine,  and  radionuclide  are 
approximated  via  standard  Galerkin  procedures.  Time  stepping  can  be 
accomplished  via  implicit  finite  difference  or  modified  method  of 
characteristics  techniques  for  the  transport-dominated  equations. 
Optimal  order  error  estimates  are  derived  when  the  imposed  external 
flows  are  smoothly  distributed.  Concepts  of  adaptive  local  grid 
refinements  are  also  discussed. 


When  things  go  wrong... 

P  W  Gaffney,  IBM,  Norway. 

This  talk  is  a  personal  account  of  my  observations  on  the  state  of 
mathematical  software  and  what  I  believe  numerical  analysts  must  do 
in  order  to  influence  computer  users  to  use  good  numerical  techniques. 


Inexact  Solutions  and  Domain  Decomposition 
G  H  Golub,  Stanford  University,  U.S.A. 

In  many  situations,  domain  decomposition  is  equivalent  to  partitioning  a 
matrix  into  submatrices  and  then  solving  a  system  of  equations  associated 
with  each  submatrix.  Unfortunately,  it  is  not  possible  to  solve  each 
subsystem  exactly.  An  inner-outer  iteration  procedure  is  necessary  in 
obtaining  the  approximate  solution  on  the  subdomain.  In  this  talk,  we 
analyse  the  Chebyshev  semi-iterative  method  for  matrices  that  possess 
"Property  A"  and  show  the  co  nee  properties  of  the  algorithm. 
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Oil  the  Dynamics  of  Certain  Time  Integrators  with  Step-size  Control 
D  F  Griffiths,  University  of  Dundee. 


When  evaluating  some  low  order  predictor-corrector  methodswith  step-size 
control  for  possible  use  in  tine  dependent  PDFS,  an  extremely  erratic 
behaviour  was  encountered  as  the  solution  approached  its  asymptotic  state. 
This  behaviour  could  not,  in  general,  be  alleviated  by  adding  constraints 
such  as  imposing  a  bound  on  the  ratio  of  consecutive  time  steps. 

In  this  talk  we  shall  investigate  the  causes  of  such  phenomena  in  very 
simple  situations  and  analyse  certain  schemes  both  from  the  point  of  view 
of  dynamical  systems  and  by  the  use  of  modified  equations. 


Some  Hew  Classes  of  Mixed  Finite  Element  Methods 
T  J  R  Hughes,  Stanford  University,  U.S.A. 


We  describe  two  general  classes  of  mixed  finite  element  methods  which 
possess  better  stability  properties  than  the  classical  Galerkin 
formulation.  The  first  class  enables  convergence  theorems  to  be 
established  while  circumventing  entirely  the  Babuska-Brezzi  condition. 
Methods  of  this  type  (termed  "CBB  methods”)  have  been  developed  which 
permit  a  wide  variety  of  convergent  interpolations,  many  of  which  are 
unstable  in  the  Galerkin  formulation.  For  example,  n  CBB  formulation 
of  the  Stokes  problem  is  described  which  is  convergent  for  all 
combinations  of  continuous  velocity  and  continuous  or  discontinuous 
pressure  interpolations.  The  second  class  of  methods  permits  satisfaction 
of  the  BabuSka-Brezzi  condition  and  a-ellipticity  condition  for  a  much 
wider  class  of  interpolations  than  the  Galerkin  method  (”SBB  methods”). 
Applications  to  problems  arising  in  Structural  mechanics  are  described. 


Order  Barriers  for  Difference  Schemes  for  Linear  and  Nonlinear 
Hyperbolic  Problems 
R  Jeltsch,  RWTil  Aachen,  Germany. 

For  the  constant  coefficient  advection  equation,  starting  from  known 
results  for  one-step  schemes,  ve  conjecture  the  existence  of  an  order 
Courant-Friedrichs-Lewy  condition  for  multistep  methods.  Partial 
results  supporting  this  conjecture  are  given.  In  the  nonlinear 
part  of  the  paper  we  concentrate  on  one-step  schemes  for  a  scalar 
hyperbolic  conservation  law.  The  order  conditions  sre  given  and  the 
order  of  the  Random-Choice  method  is  investigated.  Old  and  new  order 
barriers  due  to  different  properties  of  the  schemes  are  presented. 
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Arithmetic  for  Vector  Processors 
U  Kulisch,  University  of  Karlsruhe,  Germany. 


In  electronic  computers  the  elementary  arithmetic  operations  are  these 
days  generally  approximated  by  floating-point  operations  of  highest 
accuracy.  Computer  arithmetic,  however,  can  be  extended  so  that  the 
arithmetic  operations  in  the  linear  spaces  and  their  interval 
correspondents  which  are  most  conmonly  used  in  computation  can  be 
performed  with  highest  accuracy  on  digital  computers  also.  The 
basic  building  blocks  of  these  additional  arithmetic  operations  are 
discussed. 

In  recent  years  there  has  been  a  significant  shift  of  numerical 
computation  from  general  purpose  computers  towards  vector  and 
parallel  computers  -  socalled  supercomputers.  Besides  the 
elementary  arithmetic  operations  ♦,  -,  *,  /,  supercomputers  usually 
offer  a  number  of  compound  operations  as  additional  elementary 
operations.  Such  elementary  compound  operations  are  for  example: 
"multiply  and  add",  "accumulate"  or  "multiply  and  accumulate". 

Also  these  operations  shall  always  deliver  an  answer  of  highest 
accuracy  whatever  the  data  are.  Diverse  methods  are  discussed 
which  allow  a  fast  and  correct  computation  of  elementary  compound 
operations  as  well  as  all  basic  operations  for  extended  computer 
arithmetic.  Pipeline  techniques  are  used  to  obtain  high  speed  for 
these  operations.  In  a  scalar  product  computation  these  methods 
allow  the  execution  of  a  multiplication  and  the  corresponding  addition 
in  each  cycle  time. 

The  methods  discussed  in  the  paper  can  also  be  used  to  build  a  fast 
unit  for  advanced  arithmetic  for  micro  computers  in  VLSI-technology. 


Robust  Non-linear  Data  Fitting 
K  Madsen,  Technical  University,  Denmark. 


During  the  last  decades  there  has  been  much  interest  in  replacing  the 
least  squares  criterion  in  data  fitting  with  a  more  “robust"  criterion. 
One  of  these  is  Huber  criterion  where  the  smaller  residuals  are 
measured  by  the  l^-norm  and  the  larger  by  the  1^-norm. 

Me  give  an  introduction  to  the  Huber  criterion  for  data  fitting. 
Algorithms  for  the  linear  case  are  reviewed,  and  a  new  trust  region 
method  for  the  nonlinear  case  is  presented.  The  latter  method  is 
illustrated  through  numerical  examples. 


Subdivision  and  Itef inement  Algorithms  for  Curves  and  Surfaces 
C  Micchc lli,  IBM,  New  York. 


We  will  survey  some  subdivision  techniques  used  in  computer  aided  design 
for  the  computation  of  curves  and  surfaces  and  also  report  on  recent 
joint  work  with  li.  I’rautsch. 
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Multivariable  Approximation  Using  Radial  Basis  Functions 
M  J  D  Powell,  University  of  Cambridge. 


Our  approximating  functions  have  the  form 
a 

a(x)  -  E  A.$(||  x-x.  ||  _),  xtlRn 
i-1  1  1  z 

wliere  the  (A,  :  i*l,2,+. ..  ,ro}  are  real  coefficients,  4(.)  is  a 
fixed  function  from  1R  to  JR  ,  and  [x.  :  i**l,2, . . .  ,m)  is  a 
set  of  prescribed  points  in  ]Rn  .  They  seem  to  be  highly  useful 
for  approximation  to  continuous  functions  from  1R  to  R, 
because  the  points  (x.  :  i-l,2,...,m}  can  be  made  more  dense 
where  higher  accuracy  or  more  flexibility  is  required.  Ue 
consider  published  work  on  choices  of  f(.),  interpolation  on 
(x.  :  i-1,2, ... .m} ,  basis  functions  for  the  linear  space  that  is 
spanned  bv  {$(  ||  x-x.  ||  )  s  xelR"}  (i-1,2, ...  ,ra),  the  addition  of 
a  low  ordec  polynomial  to  s(.),  and  convergence  properties.  The 
cases  4(r)-r  and  4(r)"r^  when  n-2  are  demonstrated  by  nmerical 
examples. 


Review  of  some  Results  on  Modulations!  Instabilities  in  Nonlinear 
Difference  Schemes 

D  M  Sloan*  and  B  M  Herbst,  University  of  Strathclyde. 


The  numerical  solution  of  a  nonlinear  difference  equation  may  become 
unstable  in  situations  where  the  exact  solution  of  the  difference 
equation  should  remain  bounded.  Newell  discussed  this  phenomenon 
and  attributed  it  to  a  focusing  siechanism  which  localises  energy  on 
the  spatial  grid.  Here  we  review  some  results  obtained  using 
leap-frog  discretisations  of  the  advection  equation  and  the  KdV 
equation.  numerical  results  show  a  modulation  of  the  envelope  of 
the  basic  solution  and,  as  this  modulation  takes  place,  there  is  a 
growth  in  Fourier  modes  in  the  vicinity  of  the  fundamental  mode. 

The  solution  evolves  in  time  with  variations  which  depend  on  more 
than  one  tisie  scale.  A  multiple  scales  analysis  is  described  which 
throws  some  light  on  the  structure  of  the  instability. 


A  Study  of  Rosenbrock  Methods  for  Differential  and  Differential- 
Algebraic  Problems 

G  Wanner,  University  of  Geneva,  Switzerland. 


7 


Computational  Method a  for  Bifurcation  Problems  with  Symmetric!  and 
Applications  to  Steady  States  of  n-box  Reaction-diffusion  Models 
B  Werner,  University  of  Hamburg,  Germany.  ~ 


Ue  consider  parametrized  nonlinear  equations  (equilibrium  problems) 
(1.1)  g(x,A)  -  0,  g:  X  «  ®  +  X,  X  - 


satisfying  a  symmetry  condition 


(1.2)  g(yx,A)  -  Yg(x,A),  xeX,  1  e  S,  1  t  T, 

with  a  symmetry  group  T  c  0(n),  Discretized  nonlinear  boundary 
value  problems  with  syanetric  domain  or  symmetrically  coupled  n-box 
reaction-diffusion  models  lead  to  (1.1),  (1.2). 

Due  to  the  symmetry,  (1.1)  can  carry  a  rich  bifurcation  structure: 
there  might  be  several  E-branches  of  solutions  of  (1.1)  having  the 
symmetry  of  a  subgroup  E  of  T  and  being  connected  with  other 
E-branches  via  E-E-(syrametry-)  breaking  (multiple)  bifurcation 
points. 

It  is  shown  how  group  theoretical  ideas  and  bifurcation  results 
(Cicogna  [1981],  Vanderbauwhede  [1982),  Co lubitsky -Stewart  -  Schaeffer 
[1987)  allow  an  efficient  computation  of  the  solution  diagram  of  (1.1), 
especially  of  the  bifurcation  points.  The  essential  step  is  a 
reduction  process  corresponding  to  suitable  isotropy  subgroups  E  of 
T.  Assuming  irreducible  representations  of  T  on  the  kernel  of  a 
singular  point  of  g,  multiple  bifurcation  points  are  reduced  to 
simple  (e.g.  pitchfork)  bifurcation  points  and  the  extended  system 
in  Werner-Spence  11989]  (treating  the  case  of  Z2-symraetry)  can  be 
used  for  their  computation.  Numerical  results  for  certain  9-box 
Brusselator  models  are  given.  Related  aspects  for  the  construction 
of  Gaussian  cubature  rules  for  certain  domains  with  synsnetry  are 
mentioned. 
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SECTION  B 


Some  Practical  Experience  on  a  Parallel  Computer 

C  A  Addison*;  J  Petersen;  R  H  Chamberlain,  Chr.  Michelsen  Institute, 
Bergen,  Norway. 


He  present  an  overview  of  the  parallel  computing  activities  at  Chr. 
Michelsen  Institute  over  the  last  year  and  a  half.  These  include 
aspects  of  seismic  processing,  the  solution  of  the  2D  wave  equation 
in  inhomogeneous  media  and  the  solution  of  dense  systems  of  linear 
equations.  This  presentation  focuses  on  the  latter  two  topics. 

Certain  operations  among  the  nodes  of  the  Hypercube  are  used  frequently. 
A  high  level  library  of  such  operations  is  evolved  in  order  to  make 
using  the  Hypercube  easier.  It  is  intended  that  such  a  library  hide 
machine  dependencies  from  users,  thua  providing  a  common  interface  for 
an  application  of  any  Hypercube  based  architecture. 

He  conclude  with  a  susaaary  of  what  we  have  learned  in  these  initial 
efforts  of  solving  problems  in  parallel  and  provide  an  indication  of 
what  our  future  developments  might  be. 


Experiments  with  the  BFCS  and  DFP  Ppdatea  on  Various  Positive  Definite 
Matrices 

M  Al-Baali,  University  of  Daauscus,  Syria. 

In  this  talk  we  consider  the  possibility  of  correcting  the  current  matrices, 
which  approximate  the  Hessian  matrix  or  its  inverse,  in  certain  sense, 
before  updating  either  by  the  BFCS  or  by  the  DFP  formulae.  This  approach 
leads  to  definite  various  classes  of  formulae  for  minimizing  unconstrained 
problems,  and  one  of  them  is  the  well  known  Broyden  family  of  formulae. 
Numerical  results  are  given  to  illustrate  most  cases. 


A  Comparison  of  One-Dimensional  Solution  Techniques  Applied  to  the  2-D 
Euler  Equations 

J  J  Barley,  University  of  Reading. 


Approximate  solutions  of  the  two-dimensional  Euler  equations  of  compressible 
flow  are  often  obtained  by  operator  splitting  where  a  one-dimensional 
difference  operator  is  applied  alternately  in  each  coordinate  direction. 
Alternatively,  the  simultaneous  application  of  one  dimensional  difference 
operators  in  each  direction  nay  be  used. 

A  comparison  is  presented  between  such  techniques  for  two  test  problems 
(mild  and  severe)  using  both  the  flux-difference  split  scheme  of  Roe  and 
the  flux-vector  split  scheme  of  Steger-Varming.  Results  are  illustrated 
graphically. 
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Natural  Conti  js  Extension*  of  Runge-Kutta  Method!  for  V  ;rra 
Integral  Equations  of  the  Second  Kind  and  their  Application* 

A  Be  lien;  Z  Jackiewicsj  R  Vermiglio*;  H  Zennaro,  Univeriita'  di 
Udine,  Italy. 


Ue  consider  a  very  general  class  of  Runge-Kutta  set hods  for  the  numerical 
solution  of  Volterra  integral  equations  of  the  second  kind,  which 
includes  as. special  cases  all  the  no it  important  Methods  which  have  been 
considered  in  the  literature.  The  Bain  purpose  of  this  paper  is  to 
define  and  prove  the  existence  of  the  Natural  Continuous  Extensions 
(UCEs)  of  Runge-Kutta  methods,  i.e.  piecewise  polynomial  functions  which 
extend  the  approximation  at  the  grid  points  to  the  whole  interval  of 
integration.  The  particular  properties  required  of  the  NCEs  allow  us 
to  construct  the  tail  approximations,  which  are  quite  efficient  in  terms 
of  kernel  evaluations  and  also  result  in  good  stability  properties  of 
the  proposed  numerical  methods. 


Clobal  Error  Estimation  in  the  Method  of  Lines  for  Parabolic  Equations 
M  Berlins,  Leeds  University. 


A  method  is  described  for  obtaining  an  indication  of  the  error  in  the 
numerical  solution  of  parabolic  partial  differential  equations  using 
the  method  of  lines.  The  error  indicator  is  derived  by  using  a 
combination  of  existing  global  error  estimating  algorithms  for  initial 
value  problems  in  ordinary  differential  equations,  [1]  with  estimates 
for  the  p.d.e.  truncation  error  which  are  obtained  by  using  Richardson 
extrapolation.  An  implementation  of  the  algorithm  is  described  and 
numerical  examples  are  used  to  illustrate  the  reliability  of  the  error 
estimates  that  are  obtained.  The  efficiency  of  the  algorithm  is 
evaluated  and  various  different  approaches  compared.  The  application 
of  the  method  to  the  design  of  reliable  algorithms  for  partial 
differential  equations  is  considered. 

[1]  L.F.Shampine:  Globa)  Error  Estimation  for  Stiff  O.D.E.s.  Proc.  of 
Dundee  Conference  on  Numerical  Analysis  1983,  Springer  Verlag  Lecture 
Notes  in  Mathematics  No  1066. 


Numerical  Hethods  for  the  Space  Charge  Equation 
C  Budd,  Oxford  University. 


Die  space  charge  equation  is  a  third  order  nonlinear  PDE  of  ^xed  type 
which  models  the  steady  electrostatic  field  in  a  neighbourhood  of  a 
coronating  electrode.  By  employing  a  hodograph  transformation  we  may 
rewrite  the  equation  as  a  coupled  nonlinear  elliptic  equation  and  an 
ODE  on  a  rectangular  domain.  This  transformation  leads  to  an  effective 
numerical  scheme  based  upon  a  finite  difference  approach,  which  can  cope 
with  a  wide  variety  of  geometries  and  boundary  conditions  for  the  original 
problem. 
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Convergence  of  Univariate  Quasi-Interpolation  Using  Hultiquadrics 
M  D  Buhmann,  University  of  Cambridge. 

Univariate  interpolation  using  the  radial-basis-function  R  -*•  R 
for  an  infinite  number  of  equally-spaced  real  data-points  of  spacing 
h  >  0  can  be  performed  by  choosing  the  element  of  the  linear  space 
X  spanned  by  (♦(|.~jh|)  j  j  t  Z)  that  satisfies  f»(ih)  »  6.-  and 
defining  the  interpolant  s:  R  -*•  T  to  ft  R  ■*  R  by  11 

09 

s(x)  -  Z  f (jh)<i(x-jh),  x  e  R, 

jm—cn 

which  yields  s(ih)  »  f(ih)  V  i  e  Z. 

Restricting  f  to  a  finite-dimensional  subspace  of  X  and  giving  up 
the  condition  {i(ih)  ■  6^q,  we  consider  quasi-interpolants  s:  R  ■+•  R 
of  the  form  given  and  study  their  convergence  behaviour  with  regard  to 
localization  properties  of  ^i.  If  twice  differentiable  functions  f 
with  |f(x)|  »  0(|x|)  for  x  ±t  °°  are  approximated  and  if  Hardy's 
multiquadric  formula  4(r)  ■  /r!  ♦  c*  is  chosen  for  c  >  0,  localization 
is  sufficient  to  prove  convergence  of  the  approximation  scheme  as  h  0. 
When  fir)  •  l//r*  ♦  c  ,  however,  as  in  the  inverse  multiquddrics 
formula,  sufficient  localization  cannot  be  obtained. 

Studying  the  influence  of  the  choice  of  the  parameter  c  in  the 
multiquadrics  formula  on  the  quality  of  the  approximation,  we  find 
that  c/h  has  to  remain  bounded  as  h  +  0.  Otherwise  there  are 
smooth  and  bounded  functions  for  which  the  quasi-interpolants  given 
do  not  converge  to  the  said  functions. 

He  finally  give  numerical  examples  that  suggest  that  in  the  cases  of 
nltiquadrics  and  inverse  multiquadrics  where  now  again  <i(ih)  » 
and  f  is  an  element  of  the  whole  space  X,  good  localization 
properties  of  f  are  achieved. 


P**  Characteristic  Galerkin  Method  for  Nonlinear  Hyperbolic  Systems 
P.  Childs,  Oxford  University  Cot^uting  Laboratory. 


He  describe  the  application  of  the  characteristic  Galerkin  method  to 
nonlinear  systems  of  hyperbolic  conservation  laws,  especially  for  shock 
modelling.  For  a  scalar  equation,  the  basic  characteristic  Galerkin 
scheme,  based  on  piecewise  constant  elements  is  first  order;  adaptive 
recovery  techniques  give  second  order  accuracy  for  smooth  solutions. 

The  method  is  extended  to  systems  via  flux-splitting.  Our  formulation 
easily  allows  the  use  of  moving  and  adaptive  grids.  There  is  no 
stability  restriction  imposed  on  the  timestep:  choice  of  the  grid  and 
timestep  can  be  made  solely  from  considerations  of  accuracy.  A  full 
analysis  for  a  scalar  equation  demonstrates  several  interesting  features 
of  the  method.  We  illustrate  the  application  to  systems  of  equations 
through  a  series  of  one-dimensional  test  problems. 
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A  Numerical  S  ,  of  Nonlinear  Waves  In  a  Fluidized  Bed 
I  Christie*;  u  Ganser,  West  Virginia  University,  U.S.A. 

A  third  order  nonlinear  problem  approximating  the  vertical  flow  in  a  one 
dimensional  fluidized  bed  is  solved  by  a  Petrov-Galerkin  method. 

Piecewise  linear  trial  functions,  cubic  spline  teat  functions  and  product 
approximation  lead  to  a  fourth  order  discretization  in  space.  Second 
order  backward  differentiation  is  used  for  the  time  integration. 
Experiments  show  that  while  some  values  of  the  time  step  may  produce 
stable  results,  a  reduction  in  the  time  step  introduces  instability. 

This  behaviour  is  confirmed  by  a  von  Neumann  stability  analysis  of  a 
simplified  problem. 


Numerical  Methods  for  y1*  -  f(x,y)  via  Rational  Approximations  for  the 
Cosine 

J  P  Coleman,  University  of  Durham. 


Many  numerical  methods,  when  applied  to  the  scalar  test  equation  y"  “  — oj^y, 
yield  a  three-term  recurrence  relation  which  has  a  characteristic  equation 
of  the  form 


-  2R  (\)*K  +  1-0; 
nm 

here  v  -  (ah,  where  h  is  the  steplength,  and  Rm  is  a  rational  function. 
Examples  include  Runge-Kutta-Ny strom  methods  [van  der  Houwen  and  Somaeijer 
(1987a)].  Numerov's  method  and  several  related  hybrid  methods  "Chawla 
(1983,  1984),  Chawla  and  Ran  (1984,  1986)),  and  some  other  multistep  methods 
[Cash  (1981),  Thomas  (1984),  van  der  Houwen  and  Sotmneijer  (1987b)]. 

The  properties  of  Rnm(vl)  as  a  rational  approximant  for  cos  v  can  be  used  to 
investigate  the  order  of  dispersion  (phase-lag)  and  the  stability  of  such 
methods.  In  particular: 

(1)  For  given  n  and  m,  the  degrees  of  numerator  and  denominator,  the  maximal 
order  of  dispersion  is  2n  +  2m,  which  is  obtained  when  R^lv2)  i*  the 
[2n/2m]  Pade  approximant  for  cos  V. 

(2)  For  explicit  methods  of  maximal  order  of  dispersion  R^iv2)  -  P20(v)> 
the  Maclaurin  polynomial  of  degree  2n  for  cos  v.  Any  such  method  has 
an  interval  of  periodicity  (0,v2n)  where 

v.  <  v  <  u  <  ...  <  it  <  v  <  u  <  ...  <  2ir. 

•  is  2  0 

(3)  Methods  corresponding  Co  the  [0/2m]  Fade  approximants  for  cos  V 
are  P-s table. 
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Calculate...  of  Solutions  of  some  P.D.E.  with  Two-  ano  ..ireedimensional 
Singularities 

L  Collatz,  University  of  Hamburg,  Germany. 


Inclusion  theorems  for  solutions  of  some  types  of  boundary  value  problems 
with  singularities  in  and  R^.  Improvement  of  the  error  bounds 

with  aid  of  Bessel  functions  of  not  integer  index. 


Primal  Methods  are  Better  than  Dual  Methods  for  Solving  Overdetermined 
Linear  Systems  in  the  1-infinity  Sense- 

A  R  Conn; Tuying  Li;  R  Bartels:  University  of  Waterloo,  Canada. 

It  is  generally  believed  that  methods  based  upon  solving  the  dual  of  the 
linear  programing  formulation  of  the  linear  Chebychev  problem  are 
superior  to  methods  based  upon  the  primal  formulation. 

We  present  a  primal  approach  along  with  evidence  that  for  random  problems 
primal  methods  are  to  be  preferred.  We  then  show  that  the  method  that 
is  generally  considered  the  best  algorithm  for  discrete  linear  approximation 
in  the  1-infinity  norm  actually  owes  much  of  its  efficiency  to  effective 
choice  of  a  starting  point.  this  phenomenon  is  explained  with  the  aid 
of  results  from  classical  minimax  theory,  which  is  then  used  to  determine 
a  suitable  initial  point  for  a  primal  method. 

With  this  enhancement  it  would  appear  that  the  primal  approach  is  also 
preferable .for  data-fitting  problems. 


Implementation  Schemes  for  Implicit  Runge-Kutta  Methods 
G  J  Cooper,  University  of  Sussex. 


Iterative  procedures,  for  the  implementation  of  implicit  methods  applied 
to  stiff  problems,  are  considered.  Some  schemes,  which  are  based  on  the 
use  of  projection  methods,  are  obtained.  These  schemes  avoid  the  need 
to  solve  systems  of  linear  equations. 
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Least  Squares  .rcle  Fitting  to  Data  with  Specified  Uncert. _ Vty  Ellipses 

N  G  Goxj  H  M  Jones*,  Rational  Physical  Laboratory. 


We  consider  the  problem  of  fitting  a  circle  to  a  set  of  data  points 
specified  in  terms  of  their  Cartesian  coordinates.  We  assume  that 
the  data  adequately  represents  a  circle  and  that  associated  with  each 
point  there  is  an  uncertainty  ellipse  describing  the  measurement 
error.  We  formulate  a  weighted  nonlinear  least  squares  problem  in 
order  to  determine  unbiased  estimates  of  the  coordinates  of  the  centre 
and  the  radius  of  the  circle  which  best  fits  the  data.  The  weights  used 
in  forming  the  objective  function  for  this  optimization  are  derived  from 
the  given  ellipse  parameters.  The  resulting  problem  displays  structure 
which  is  exploited  when  the  Causs-Newton  algorithm  is  used  to  obtain  a 
solution.  Initial  values  for  the  method  are  estimated  from  a  linear 
approximation  to  the  nonlinear  fitting  problem. 

In  addition  to  estimates  of  the  circle  parameters  the  algorithm  produces 
their  variance /covariance  matrix  from  which  an  error  ellipse  for  the 
centre  of  the  circle  may  be  deduced. 

Problems  involving  circle  fitting  have  arisen  recently  in  two  distinct 
areas  of  work  at  NPL  -  dimensional  metrology  and  microwave  measureaunt. 

A  data  set  from  one  of  these  areas  is  used  to  illustrate  the  algorithm. 


Stability  of  the  Exact  Collocation  Method  for  the  Second  Kind  Volterra 
Integral  Equation  with  Degenerate  Kernel 
M  R  Crisci;  E  Russo;  A  Vecchio,  University  of  Naples. 

Ue  analyze  the  stability  properties  of  the  exact  collocation  method 
applied  to  a  second  kind  Volterra  integral  equation  with  degenerate 
kernel: 

t 

(1)  y(t)  -  g(t)  ♦  (  Z  A  <t)  B, (t)y(s)ds 

I  i-1  1 

Jo 

The  kernel  of  (1)  is  a  generalization  of  polynomial  convolution  kernel, 
and  it  can  be  considered  a  significant  test  ternel,  as  the  class  of  (1) 
is  dense  in  the  class  of  all  continuous  kernels  k(t,s). 

We  derive  a  recurrence  formula  for  a  vector  containing  the  numerical 
solution  in  the  collocation  points,  and  we  obtain  stability  conditions. 

Moreover,  for  the  known  convolution  test  equation  whose  kernel 

(2)  k(t,s)  ■  \  *  U  (t-s) 


is  a  particular  case  of  (l),  we  individualize,  under  the  hypothesis 
that  the  collocation  points  are  symmetric,  the  stability  regions. 
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Automatic  Imple.  .ation  of  an  Initial  Value  Method  for  use  >wo-Point 
Boundary  Problems  ” 

R  England*;  R  M  Mattheij,  Polytechnic  of  the  South  Bank. 


Two-point  boundary  value  problems  often  have  sharp  boundary  layers,  or 
rapidly  varying  fundamental  modes.  These  make  a  uniform  discretization 
inefficient,  and  iterative  adaptive  refinement  may  also  be  undesirable, 
particularly  for  linear  problems.  With  a  sequential  marching  approach, 
it  should  be  possible  to  select  reasonable  step  sizes  in  a  single  pass, 
based  on  numerical  estimation  of  layers  and  of  a  smooth  particular  solution. 
Normal  integrators  for  initial  value  problems  select  their  step  sizes  on 
the  basis  of  asymptotic  error  estimates  for  small  step  sizes.  However, 
these  may  be  large  in  the  presence  of  fast  growing  fundamental  modes, 
unless  the  step  size  is  much  smaller  than  necessary.  In  this  talk,  a 
new  form  of  error  indicator  will  be  outlined,  which  enables  an  integrator 
to  take  an  appreciable  number  of  steps  of  a  suitable  size  for  such  a 
problem,  without  generating  excessively  unstable  solutions.  After  a 
starting-up  stage,  a  special  explicit  predictor  may  be  matched  to  an 
appropriate  implicit  corrector  in  such  a  way  that  their  difference,  while 
asymptotically  correct  for  small  step  sizes,  is  also  an  appropriate  error 
indicator  for  large  step  sizes.  Efficient  and  stable  iterative  methods 
will  also  be  discussed  for  solving  the  linked  correction  formulae  which 
are  used,  while  maintaining  their  stability  properties.  A  few  results 
will  show  the  feasibility  of  finding  smooth  solutions  with  appropriate 
step  sizes. 


An  Analysis  of  Defect-Control  Strategies  for  Runge-Kutta  Methods 
W  II  Enright,  University  of  Toronto,  Canada. 

Over  the  last  three  decades  there  has  been  a  considerable  change  in 
what  is  expected  of  a  general  purpose  Runge-Kutta  method  for  initial 
value  problems.  Discrete,  f ixed-stepsize  methods  have  been  replaced 
by  discrete,  variable-step  methods  which  attempt,  with  mixed  success, 
to  match  the  achieved  accuracy  with  a  tolerance  parameter,  TOL.  Now, 
with  the  development  of  inexpensive  interpolants,  it  is  feasible  to 
expect  a  Runge-Kutta  method  to  produce  a  continuous  solution 
approximation  with  an  accuracy  which  is  directly  related  to  TOL.  In 
this  paper  we  analyse  several  error  control-strategies  which  are 
designed  to  achieve  this  objective.  We  will  present  theoretical  and 
experimental  evidence  to  show  that  this  approach  is  viable  and  can  lead 
to  a  family  of  methods  that  are  effective  and  much  easier  to  use  than 
existing  methods. 
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Pi.8crete-ti.ae  Finite  Element  Methods  for  Parabolic  Integto-Differential 
Equations 

G  Fair-weather*;  E  G  Yanik,  University  o{  Kentucky,  U.S.A. 

In  the  modeling  of  some  space-time  dependent  physical  systems  in  such 
fields  as  heat  transfer,  nuclear  reactor  dynamics  and  viscoelasticity, 
it  is  often  necessary  to  take  into  account  memory  effects.  In  the 
mathematical  model  of  such  a  system,  this  results  in  the  inclusion  of 
an  integral  term  in  the  basic  differential  equation  yielding  a  partial- 
integro  differential  equation  (FIDE).  The  equations  considered  in  this 
paper  are  PIPES  of  parabolic  type  of  the  form 

t 

(l)  u£  ♦  V(a(x,t,u)Vu)  »  jf (x,t,s,u(x,s),Vu(x,s))ds,  (x, t)  t $J  x (0,T], 

0 

where  Q  is  a  bounded  domain  in  R,  n-1,2.  While  much  attention  has 
been  paid  in  the  recent  literature  to  the  study  of  various  properties  of 
solutions  of  problems  of  this  kind,  little  has  been  reported  on  the 
development  and  analysis  of  numerical  techniques.  In  this  paper,  we 
describe  several  discrete-time  finite  element  procedures  for  the  numerical 
solution  of  (1).  Galerkin's  method  or  spline  collocation  at  Gauss  points 
is  used  for  the  spatial  discretization,  and  both  nonlinear  and  algebraically 
linear  time-stepping  methods  are  discussed.  Under  certain  conditions,  the 
methods  yield  approximations  which  are  second  order  accurate  in  time  and  of 
optimal  accuracy  in  space.  Similar  methods  for  certain  hyperbolic  PIDEs 
are  also  described. 


An  Error  Bound  for  Multivariate  Lacunary  Interpolation 
T  Fawzy,  Suez-Canal  University,  Egypt. 


In  Cl],  we  presented  a  method  for  solving  the  general  lacunary  interpolation 
problem,  using  polynomial  splines,  which  gives  optimal  error  bounds.  In 
[2],  we  introduce  a  tensor  version  of  our  general  method  which  can  be  used 
to  solve  lacunary  interpolation  problems  in  two  or  more  variables. 

The  purpose  of  this  paper  is  to  obtain  an  optimal  error  bounds  for  a  special 
case  of  the  method  described  in  [21. 

Cl]  Th.  Fawzy  and  L.L.  Schumaker:  A  piecewise  polynomial  lacunary 
interpolation  method.  Journal  of  Approximation  Theory,  Vol.48, 

Ho. A,  December  1986,  pp. 407-426. 

[2]  Th.  Fawzy  and  L.L. Schumaker :  Multivariate  Lacunary  Interpolation. 
Proceeding  of  the  Fifth  International  Symposium  on  Approximation 
Theory  held  at  Texas  A  4  M  University  on  January  13-17,  1986. 

Academic  Press  1986,  pp. 339-342. 
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Nonlinear  Programming  an<f  Nonsmooth  Optimization  by  Successive  Linear 
Programming 

R  Fletcher*;  E  Sainz  de  la  Maza,  University  of  Dundee. 

Methods  are  considered  for  solving  nonlinear  programing  problems  using 
an  exact  penalty  function.  LP-like  subproblems  incorporating  a 

trust  region  constraint  are  solved  successively  both  to  estimate  the 
active  set  and  to  provide  a  foundation  for  proving  global  convergence. 
In  one  particular  method,  second  order  information  is  represented  by 
approximating  the  reduced  Hessian  matrix,  and  Coleman-Conn  steps  are 
taken.  A  criterion  for  accepting  these  steps  is  given  which  enables 
the  superlinear  convergence  properties  of  the  Coleman-Conn  method  to 
be  retained  whilst  preserving  global  convergence  and  avoiding  the 
Maratos  effect.  The  methods  generalize  to  solve  a  wide  range  of 
composite  nonsmooth  optimization  problems  and  the  theory  is  presented 
in  this  general  setting.  A  range  of  numerical  experiments  on  small 
test  problems  is  described. 


The  Pgeudosgectral  Method:  Comparisons  with  Finite  Differences 
for  the  Elastic  Wave  Equation 
B  Fornberg,  Exxon  Research,  U.S.A. 


The  pseudospectral  method  is  a  technique  for  solving  hyperbolic  systems 

of  equations.  It  is  much  leas  flexible  than  finite  differences,  for 

example  with  regard  to  boundaries  and  adaptive  meshes.  However,  in 

applications  where  its  limitations  are  not  inq>ortant,  very  large  savings 

in  computer  time  and  storage  may  be  achieved.  Early  applications  ( 

include  homogeneous  turbulence,  nonlinear  waves  and  weather/climate  ' 

modeling.  In  this  presentation,  we  will  first  briefly  review  the 

background  of  the  method  and  then  illustrate  its  performance  in  the 

case  of  elastic  waves  propagating  through  two-dimensional  discontinuous 

media. 


Adaptive  Grid  Methods  for  Nonlinear  Dispersive  Wave  Equations 
E  Fraga,  University  of  Dundee. 


Adaptive  grid  methods  for  nonlinear  dispersive  wave  equations  are  discussed. 
Particular  attention  is  placed  on  the  Korteweg-de  Vries  and  Nonlinear 
Schroedinger  equations. 
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Pictures  of  Karmarkar's  Linear  Programing  Algor ?,thm 
D  H  Gay,  AT  &  T  Bell  Laboratories,  U.S. A. 


Karmarkar's  linear  programing  algorithm  handles  inequality  constraints 
by  changing  variables  to  make  all  constraints  about  equally  distant; 
it  moves  in  the  steepest-descent  direction  seen  by  the  new  variables. 
This  paper  summarizes  four  variants  of  Karmarkar's  linear  programming 
algorithm  (primal  affine,  primal  projective,  dual  affine,  and  dual 
projective),  discusses  depicting  poly topes  (feasible  regions),  and 
presents  pictures  illustrating  the  latter  three  variants.  These 
pictures  give  an  algorithm's  eye  view  of  the  variable  changes  and 
provide  visual  verification  of  certain  properties. 


On  Approximate  Solutions  of  One  Class  of  Integral  Equations 
L  Hacla,  Technical  University  of  Poznan,  Poland. 


Various  methods  of  approxisiate  solution  for  certain  class  of  integral 
equations  are  presented.  Sufficient  conditions  for  the  convergence  of 
the  considered  methods  are  given  and  the  error  estimates  are  established. 

The  comparison  is  made  among  introduced  methods,  that  are  illustrated 
by  exanq>les. 


A  Multilevel  Parallel  Solver  for  Banded  Linear  Systems 
I  N  Ha j j ;  S  Skelboe*,  University  of  Copenhagen,  Denmark. 


The  talk  describes  an  algorithm  for  the  parallel  LU-factorization  and 
solution  of  banded  systems.  For  a  system  of  dimension  n  with  bandwidth 
2m+l,  a  parallel  computer  with  n/(2m)  processors  and  suitable  communication 
topology  will  be  able  to  perform  factorization  and  solution  in 
0(m3  log2(n/m)  )  time  units.  The  hypercube  topology  is  an  example  of  a 
suitable  communication  topology.  The  speedup  cf  the  algorithm  as  compared 
to  the  direct  solution  using  a  single  serial  processor  is  O((n/m)/log2(n/m)). 

The  algorithm  is  based  on  a  reordering  strategy  also  known  from  nested 
dissection  where  a  set  of  separator  variables  and  equations  are  identified 
and  ordered  last  to  make  the  system  appear  as  a  set  of  connected  sub-systems. 
The  sub-systems  are  LU-factored  in  parallel  leaving  the  separator  equations 
partially  factored  with  the  structure  of  a  block  tridiagonal  system.  The 
dissection,  reordering  and  parallel  partial  LU-factorization  is  applied  to 
the  block  tridiagonal  system  and  repeated  until  the  system  is  completely 
factored.  This  part  is  closely  related  to  cyclic  reduction  of  a  block 
tridiagonal  matrix.  However,  the  result  of  the  proposed  algorithm  is 
identical  to  a  standard  LU-factorization  applied  after  the  reorderings. 

The  solution  phase  can  also  be  done  in  parallel  either  after  the  completion 
of  the  LU-factorization  or  with  the  forward  substitution  pipelined  immediately 
after  the  LU-factorization. 

The  proposed  algorithm  does  not  permit  pivoting  for  numerical  reasons  but  is 
otherwise  rather  flexible.  It  can  be  applied  to  any  matrix  with  a  block 
tridiagonal  envelope  which  includes  band  matrices,  and  if  the  matrix  is 

sn<i  positive  definite,  Cholesky's  method  can  readily  be  substituted 
for  the  LU-factorization. 


Numerical  Simulation  of  Flow  in  Porous  Media 
D  E  Heath,  University  of  Reading. 


A  production  method  for  the  extraction  of  heavy  oil  is  to  inject  steam 
into  the  reservoir  so  as  to  lower  the  viscosity  of  the  oil,  thus  allowing 
extraction  by  conventional  means,  the  tracking  of  the  steam  front  during 
the  injection  cycle  being  an  important  factor  in  optimizing  the  production. 

In  this  talk  a  mathematical  model  for  the  sisqiler  problem  of  injecting 
hot  water  into  a  porous  medium  is  described  and  the  resulting  differential 
equations  are  solved  by  fiaite  element  methods. 


On  the  Growth  of  Simplices  in  Nelder-Head  Simplex  Reflections 
D  Hensley;  P  W  Smith*;  D  J  Wood3,  IMSL,  U.S.A. 

We  derive  relationships  between  the  elements  of  a  sequence  of  simplices, 
where  each  simplex  in  the  sequence  is  produced  by  applying  a  simplex 
reflection  operation  to  the  previous  element  of  the  sequence.  The 
reflection  operation  is  one  of  four  simplex  operations  used  in  the 
tlelder-Mead  simplex  algorithm  for  function  minimization.  The  relationships 
derived  here  are  of  use  in  proving  convergence  properties  of  the  Nelder-Mead 
algorithm. 


On  the  Stability  of  nonlinear  Difference  Schemes 
B  M  Herhsc,  University  of  the  Orange  Free  State,  South  Africa. 

This  calk  looks  into  more  detail  at  one  of  the  situations  described  by 
D.M.  Sloan.  It  is  shown  how  nonlinear  difference  schemes  may  be 
destabilised  through  resonances  between  the  different  waves  constituating 
the  solution.  The  multiple  scales  analysis  is  described  and  numerical 
examples  are  given 


Analysis  of  the  Enright-Kamel  Partitioning  Method  for  Stiff  ODEs 
D  J  lligham.  University  of  Manchester. 


The  use  of  implicit  formulae  in  the  solution  of  stiff  ODEs  gives  rise  to 
systems  of  nonlinear  equations  which  are  usually  solved  iteratively  by  a 
modified  Hewton  scheme.  The  linear  alge6ra  costs  associated  with  such 
schemes  may  form  a  substantial  part  of  the  overall  cost  of  the  solution. 
Enright  and  Kamel  [7,  10]  attenq)t  to  reduce  the  cost  of  the  iteration  by 
automatically  transforming  and  partitioning  the  system.  We  provide  new 
theoretical  justification  for  this  method  in  the  case  where  the  stiff 
eigenvalues  of  the  Jacobian  matrix  used  in  the  modified  Newton  iteration 
are  small  in  ntsnber  and  well  separated  from  the  other  eigenvalues.  The 
theory  of  Saad  [15,  16]  is  introduced  and  adapted  to  show  that  the  method 
uses  the  projection  of  the  Jacobian  onto  a  Krylov  subspace  which  virtually 
contains  the  dominant  subspace.  This  is  shown  to  have  favourable 
consequences.  Numerical  evidence  is  provided  to  support  the  theory. 
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Fortran  Codes  for  Estimating  the  One-norm  of  a  Real  or  Complex  Hatrix, 
with  Applications  to  Condition  Estimation 
N  J  Higham,  University  of  Manchester. 


Fortran  77  codes  SONEST  and  CONEST  are  presented  for  estioating  the  1-norm 
of  a  real  or  comp lex  matrix  respectively.  The  codes  are  of  wide 
applicability  in  condition  estimation  since  explicit  access  to  the  matrix, 
A,  is  not  required;  instead,  matrix-vector  products  Ax  and  ATx  are 
computed  by  the  calling  program  via  a  reverse  comunication  interface. 

The  algorithms  are  based  on  a  convex  optimisation  method  for  estimating 
the  1-norm  of  a  real  matrix  devised  by  Hager  [SIAM  J.  Sci.  Stat.  Coraput. 

5  (1984),  pp. 311-316].  He  derive  new  results  concerning  the  behaviour 

of  Hager's  method,  extend  it  to  complex  matrices,  and  make  algorithmic 
modifications  in  order  to  improve  the  reliability  and  efficiency. 


Numerical  Fixed  Dosain  Mapping  Solution  of  Free  Surface  Flows  Coupled 
with  an  Evolving  Interior  Field 

RCA  Hindmarsh*;  K  Hutter,  Edinburgh  University. 

The  problem  of  a  gravity-driven  free-surface  flow  coupled  with  the  evolution 
of  an  interior  field  is  solved  using  a  finite  difference  discretization  of  a 
supping  of  the  problem  onto  the  unit  square.  This  takes  advantage  of  the 
generally  similar  nature  of  the  geoaietry  and  internal  fields  over  a  wide 
variety  of  domain  sizes  so  as  to  avoid  excessively  fine  or  coarse 
discretizations  at  maximum  and  minimum  domain  sizes. 

The  cases  of  parabolic  and  hyperbolic  evolution  equations  for  the  internal 
domain  are  considered.  The  mapping  shows  that  the  true  nature  of  parabolic 
equations  is  usually  mixed  parabolic/elliptic. 

The  evolution  of  the  coupled  system  is  solved  by  an  implicit  marching 
scheme.  The  discretizations  in  space  and  in  time  are  second-order  accurate. 
Multi-point  upwinding  is  used  to  avoid  an  instability  arising  when  advective 
terms  are  large. 

The  evolution  equations  are  non-linear,  and  are  solved  using  a  nested 
Newton-Raphson  procedure.  The  nesting  is  achieved  by  using  successively 
better  approximations  to  the  true  evolution  equations.  The  matrix  that 
arises  is  solved  using  an  incomplete  Choleaky  factorization  as  a 
pre-conditioning  method  with  a  conjugate-gradient-like  (ORTHOMIN)  iteration 
procedure. 
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Improved  Accuracy  in  Generating  Forward  Difference  Estimates  for  the 
Solution  of  the  Optimal  Control  Problem 
M  K  Horn,  Messerschmitt-Bolkow-Blohm,  F.R.  Germany. 


In  numerical  solutions  of  the  optimal  control  problem,  the  use  of  forward 
difference  approximations  for  estimating  the  partial  derivatives  of  the 
cost  function  and  the  boundary  conditions  can  reduce  user  effort  and 
set-up  time.  These  advantages,  however,  can  be  off  set  by  increased 
computing  time  if  the  inaccuracies  arising  from  the  perturbation  size 
slow  down  the  convergence  rate.  By  restructuring  the  right  hand  aides 
of  the  ordinary  differential  equations  (ODEs),  one  can  essentially  integrate 
the  difference  between  the  perturbed  solution  and  the  nominal  solution 
to  the  same  accuracy  as  that  of  the  nominal  solution.  This  permits 
the  use  of  smaller  perturbations,  and  hence,  yields  higher  accuracies  in 
the  forward  difference  estimates.  The  "restructuring"  of  the  ODEs 
is  simplified  by  the  use  of  a  software  package  for  the  basic  functions, 
so  that  little  additional  user  effort  is  required  for  modeling  the 
equations.  Theory  and  results  are  presented  using  Runge-Kutta  formulas 
of  paired  orders  (4)5  and  7(8).  Particular  attention  is  paid  to  the 
treatment  of  time-related  design  parameters,  e.g.  initialization  of 
braking  features,  run-up  to  full  thrust,  etc.  Finally,  it  is  shown  that 
the  advantage  of  such  an  approach  may  be  finessed  by  restructuring  the 
integration  package  itself,  depending  upon  the  accuracy  desired  and  the 
word  length  of  the  computer. 


Improving  the  Stability  of  Predictor-Corrector  Methods  by  Residue  Smoothing 
V  J  van  der  Houwen*;  B  P  Sonmeijer,  Centre  for  Mathematics  and  Computer 
Science,  The  Netherlands. 


Residue  smoothing  is  usually  applied  in  order  to  accelerate  the  convergence 
of  iteration  processes.  Here,  we  show  that  residue  smoothing  can  also  be 
used  in  order  to  increase  the  stability  region  of  predictor-corrector 
methods.  We  shall  concentrate  on  increasing  the  real  stability  boundary. 
The  iteration  parameters  and  the  smoothing  operators  are  chosen  such  that 
the  stability  boundary  becomes  as  large  as  c(m,q)m2Vl  where  a  is  the 
number  of  right-hand  side  evaluations  per  step,  q  the  number  of  smoothing 
operations  applied  to  each  right-hand  side  evaluation,  and  c(m,q)  a 
slowly  varying  function  of  m  and  q,  of  magnitude  1.3  in  a  typical 
case.  Numerical  results  show  that,  for  a  variety  of  linear  and  nonlinear 
parabolic  equations  in  one  and  two  spatial  dimensions,  these  smoothed 

predictor-corrector  methods  are  at  least  competitive  with  conventional 
implicit  methods. 
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Nonlineai,  .ability  of  O.U.E.  Solveta 
A  Iserles,  University  of  Cambridge. 


The  standard  Dahlquist  nonlinear  stability  theory  for  numerical  ODE's 
aims  at  extending  the  scope  of  the  linear  theory.  This  has  led  to 
great  many  important  insights.  However,  in  this  talk  we  approach  the 
underlying  problem  from  the  opposite  point  of  view,  elucidating 
phenomena  that  are  genuinely  nonlinear  and  that  cannot  be  explained  by 
the  linear  theory. 

Our  analysis  focuses  on  the  Riccati  equation,  both  scalar  and  vector. 

He  investigate  the  maintenance  of  the  correct  asymptotic  behaviour  by 
multistep  and  Runge-Kutta  methods  and  examine  the  influence  of  various 
nonlinear  algebraic  solvers.  This  is  being  done  with  techniques  from 
the  theory  of  dynamical  systems. 

Some  of  our  results  will  be,  to  say  the  least,  surprising  to  the 
linearly  minded  -  e.g.,  it  is  strongly  advantageous  to  use  an  even 
number  of  corrector  iterations  in  the  Predictor-Corrector  implementation 
of  Forward  Euler  and  the  Trapezoidal  Rule. 


Order  of  Convergence  of  Quasi-Interpolation  UsinR  Radial  Basis  Functions 
I  R  H  Jackson,  University  of  Cambridge. 


In  this  paper  only  quasi-interpolation  on  a  regular  grid  in  Rn  will 
be  considered.  For  such  a  method  it  is  necessary  to  have  a  function 
*(x)  defined  on  Rn  for  which 


(|*(x)|dx<  “  and 
Rn 

In  this  case  the  quasi-interpolant  to 
a  regular  grid  of  spacing  h,  is 

ah(x)  -  E  f  (y)*(h  1(x-y)) 


*(x)dx  -  1. 

JRn 

a  function  f  define  on  Rn,  on 


where  the  sum  is  taken  over  all  of  the  mesh  points  of  a  regular  grid 
of  spacing  h,  one  of  which  is  at  the  origin. 

If  <p(.)  has  the  "radial  basis  function"  form 

*00  "  E  Mi  (|  x  -  x.  ||  , 
i-1  1  1 


which  is  achievable  when  n  is  odd  (Jackson,  DAMTP  1986/NA6),  it  is 
shown  that  further  conditions  may  be  put  on  {p^}  and  (x^)  that 
ensure  that,  on  an  integer  grid  in  the  quasi-interpolant  of  any 

polynomial  of  degree  at  most  2n  +  1  is  the  polynomial  itself.  From 
this  result  it  is  deduced  that,  if  f  is  sufficiently  smooth,  the 
quasi-interpolant  a^fx)  satisfies 

|ah(x)  -  f  (x)  |  -  0(li2nf2)  for  all  x  «  R2n<l. 


This  provides  a  surprising  extension  to  the  fact  that  linear  interpolation 
in  one  dimension  on  a  regular  grid  is  of  second  order. 
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EPUCOLi  A  Hort.  .f icient  PDECOL  Code  /' 

P  Keast*;  P  Muir,  Dalhuusie  University,  Canada. 

the  algorithm  PDECOL  (Algorithm  540,  Trans.  Math.  Soft.,  vol.  5, 
pp. 326-351)  is  a  popular  code  among  acientista  who  have  to  solve 
systems  of  non-linear  partial  differential  equations.  The  code 
is  based  on  a  method  of  lines  approach,  with  collocation  in  the 
space  variables  to  reduce  the  problem  to  a  system  of  ordinary 
differential  equations.  There  are  three  principal  components  : 
the  basis  functions  employed  in  the  collocation;  the  method  used 
to  solve  the  system  of  ordinary  differential  equations;  and  the 
linear  equation  solver  which  handles  the  linear  algebra.  This 
talk  will  concentrate  on  the  third  component,  and  will  report  on 
improvements  in  the  performance  of  PDECOL  resulting  from  replacing  the 
linear  algebra  modules  of  the  code  by  modules  which  take  full  advantage 
of  the  special  structure  of  the  equations  which  arise.  Savings  of  up 
to  50Z  in  total  execution  time  can  be  realised. 


Calerkin/Runge-Kutta  Discretizations  for  Parabolic  Equations 
S  L  Keeling,  NASA  Langley  Research  Centre,  U.S.A. 


A  new  class  of  fully  discrete  Galerkin/Runge-Kutta  methods  is  constructed 
and  analyzed  for  parabolic  equations  which  are  (1)  semi linear, 

(2)  linear  with  time  dependent  coefficients,  and  (3)  quasilinear. 

Unlike  any  classical  counterpart,  this  class  offers  arbitrarily  high, 
optimal  order  convergence.  In  support  of  this  claim,  error  estimates 
are  proved,  and  computational  results  are  presented.  Additionally, 
for  the  problems  in  which  the  time  stepping  equations  involve  coefficient 
matrices  changing  at  each  time  step,  a  preconditioned  iterative  technique 
is  used  to  solve  the  linear  systems  only  approximately.  Nevertheless, 
the  resulting  algorithm  is  shown  to  preserve  the  optimal  order  convergence 
rate  while  using  only  the  order  of  work  required  by  the  base  scheme  applied 
to  a  linear  parabolic  problem  with  time  independent  coefficients. 
Furthermore,  it  is  noted  that  special  Runge-Kutta  methods  allow 
computations  to  be  performed  in  parallel  so  that  the  final  execution 
time  can  be  reduced  to  that  of  a  low  order  method. 


The  Inflection  Points  on  the  Curve  of  Cubic-Quartic  Spline  Fit 
Ha-Jine  Kimn,  Ajou  University,  South  Korea. 

Since  the  slope  of  a  curve  is  maximum  and  minimum  at  the  inflection  point, 
it  is  very  interesting  to  localize  exactly  the  inflection  point  in  order 
to  know  well  the  critical  situation  of  the  given  phenomenon.  Some  of 
these  phenomena  are  known  frequently  by  the  erroneous  data  with  uncertainty. 
We  have  already  proposed  the  adoption  of  the  cubic-quartic  spline  fit  to 
represent  these  phenomena. 

We  seek  an  algorithm  for  the  systematic  usage  to  obtain  the  inflection 
points  on  the  curve  constructed  by  the  cubic-quartic  spline  fit. 
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A  Parallel,  -.lgorithn  for  the  General  LU-Factorization 
0  R  Kincaid*;  T  C  Oppe,  The  University  of  Texas  at  Austin,  U.S.A. 

An  algorithm  for  computing  in  parallel  the  general  LU-factorization  of  a 
matrix  is  presented.  As  special  cases,  one  obtains  the  Doolittle,  Crout, 
and  Cholesky  methods.  The  algorithm  was  implemented  and  tested  on  the 
Cray  X-tff/48. 


Expansion  Coefficients  for  Powers  of  gome  Trigonometric  Functions  in 
Series  of  Jacobian  Elliptic  Functions 

A  Kiper*;  S  Wrigge,  Middle  East  Technical  University,  Turkey. 

The  Fourier  series  coefficients  for  the  Jacobian  elliptic  functions  and 
their  powers  had  been  studied  and  recurrence  formulae  had  been  obtained. 
It  is  also  possible  to  obtain  expansions  for  trigonometric  functions  in 
Jacobian  elliptic  functions.  In  this  work,  the  expansion  coefficients 
for  sinB(trx)  and  cosB(irx)  in  powers  of  Jacobian  elliptic  functions, 
with  m  e  1  are  atudied.  Four-terra  recurrence  formula  are  obtained  and 
computational  aspacts  of  these  coefficients  are  discussed. 


Integral  Equation  System  for  Turbulent  Diffusion 
P  Kirkegaard*; L  Kristensen,  Riso  National  Laboratory,  Denmark. 

When  describing  the  development  of  a  Gaussian  cloud  of  particles  with 
time  in  a  turbulent  medium,  a  non-linear  two-dimensional  integral 
equation  system  of  the  Volterra  type  arises.  This  work  discusses 
an  iterative  solution  method  for  these  equations.  Semianalytical 
techniques  are  used  for  extraction  of  both  the  near-field  and  the 
far-field  properties  of  the  solution.  Coraputat ' onal  results  with 
the  algorithm  including  its  convergence  properties  are  reported. 


The  order  of  B-convergence  for  algebraically  stable  Runge-Kutta  methods 
J  F  1  M  Kraai jevanger^  University  of  Leiden,  The  Netherlands. 

In  [1]  and  [2]  the  concept  of  algebraic  stability  was  introduced  for 
Runge-Kutta  methods.  It  is  well  known  that  algebraic  stability  is 
equivalent  to  contractivity  on  the  class  of  dissipative  problems. 

In  this  talk  we  show  that  algebraic  stability  is  also  equivalent  to 
B-convergence  on  certain  classes  of  strictly  dissipative  problems. 

It  turns  out  that  the  order  of  B-convergence  is  generally  equal  to  the 
stage  order. 

[1]  K.  Burrage  and  J.C.  Butcher:  Stability  criteria  for  implicit 
Runge-Kutta  methods.  SIAM  J.  Numer.  Anal.  U>,  46-57  (19/9). 

12]  M.  Crotizeix:  Sur  la  B-stabilit£  des  methodes  de  Runge-Kutta, 

Numer.  Math.  32,  75-82  (1979). 

[3]  M.N.  Spijker:  The  relevance  of  algebraic  stability  in  implicit 
Runge-Kutta  methods.  IN:  Numerical  treatment  of  differential 
equations,  K.  Strehmel  (ed.),  Teubner  1986,  158-164, 
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On  the  Solvabil  of  the  System  of  Equations  Arising  in  Imp  Jt 
Runge-Kutta  Methods 

M  Z  Liu,  University  of  Leiden,  The  Netherlands. 


In  this  talk  we  present  a  new  condition  under  which  the  system  of 
equations  arising  in  the  application  of  an  implicit  Kunge-Kutta  method  to 
a  stiff  initial  value  problem,  have  unique  solutions.  We  show  that  our 
condition  is  weaker  than  related  conditions  presented  in  tl]  and  [21. 

It  is  proved  that  the  Lobatto  II1C  methods  fulfil  the  new  condition. 

Til  M.Crouzeix  and  W.H.  Ilundsdorfer  and  M.N.  Spijker:  On  the  existence 
of  solutions  to  the  algebraic  equations  in  implicit  Runge-Kutta 
methods.  BIT  23,  84-91  (1983). 

L 2 J  W.ll.  Huudsdorfer  and  M.N.  Spijker:  On  the  algebraic  equations  in 
implicit  Runge-Kutta  methods.  Report  NM-R  8413,  Centre  for 
Mathematics  and  Computer  Science,  Amsterdam  (1984).  To  appear 
in  SIAM  J.  Numer.Anal. 


A  Nonlinear  Multigrid  Algorithm  and  Boundary-Fitted  Coordinates  for  the 
Solution  of  a  Two-Dimensional  Flow  in  a  Branching  Channel 
C  Lonsdale*;  J  S  Bramley;  0  M  Sloan,  University  of  Bradford. 

A  recent  paper  by  Bramley  and  Sloan  (Computers  and  Fluids  (1987) 
describes  a  numerical  solution  for  two-dimensional  flow  of  a  viscous, 
incompressible  fluid  in  a  branching  channel.  A  grid  generation 
algorithm  was  used  to  map  the  solution  region  onto  a  rectangle  and  an 
upwind  finite  difference  scheme  was  then  used  to  solve  the  Navier-Stokes 
equations  in  terms  of  stream  function  and  vqrticity.  Here  we  extend 
the  earlier  work  by  describing  an  efficient  solution  of  the  problem 
using  a  nonlinear  multigrid  algorithm.  Of  particular  interest  is  the 
treatment  of  the  boundary  conditions  in  a  manner  which  does  not  destroy 
the  interior  smoothness  in  the  neighbourhood  of  the  boundary. 


A  Block-by-Block  Method  for  Volterta  Delay  InteRro-Pif ferential 
Equations  with  Weakly-Singular  Kernel  Combined  with  Product  Integration 
on  Nonuniform  Meshes 

A  Makroglou,  Ministry  of  Agriculture,  Greece. 


The  solution  of  weakly-singular  Volterra  integro-dif ferential  equations 
has  in  general  end-point  singularities  at  the  left  end-point  of  the 
interval  of  integration.  Block-by-block  implicit  Runge-Kutta  methods 
applied  to  such  equations  on  uniform  meshes  give  poor  convergence  rates. 
In  this  paper  these  methods  are  generalized  so  that  they  can  also  be 
used  with  nonuniform  mesh  and  are  applied  to  weakly  singular  Volterra 
delay  integro-diffcrential  equations  with  constant  and/or  variable  state 
independent  delay.  The  kernel  is  assumed  nonlinear  and  thus  Gauss 
integration  methods  for  the  evaluation  of  the  resulting  integrals. 

'Die  paper  is  a  generalization  of  Hakcoglou  (8,10,11,12). 
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Localise^  and  Chaotic  Patterns  in  a  two  phase  flow 
V  S  Manoranjan,  University  of  Surrey. 

A  counter-current  two  phase  flow  of  ouch  interest  is  gas-liquid  annular 
flow  in  vertical  cylinders,  A  mathematical  model  which  captures  the 
essential  features  of  this  physical  problem  is  studied  in  order  to 
understand  the  various  localised  and  chaotic  patterns  observed  in  the 
'flooding*  region. 


Improving  Gauss-Seidel  Iterations  by  Means  of  Partial  Elimination 
J  P  Milaszewicz,  Ciudad  Universitaria,  Argentina. 

If  we  consider  the  linear  system  Ax*b,  with  A  an  irreducible 
M-matrix,  the  elimination  of  some  unknowns  yields  better  convergence 
rate  for  the  corresponding  Jacobi  iterations,  when  compared  with  those 
for  the  original  system.  This  is  not  always  the  case  if  Gauss-Seidel 

iterations  are  compared.  We  give  sufficient  conditions  on  the  nodes 

to  be  eliminated  that  yield  improvement  on  Gauss-Seidel  iterations. 
Such  conditions  are  always  satisfied  when  A  is  symmetric. 


Periodic  Solutions  of  Hon  Linear  Partial  Difference  Equations  in  Reaction- 
diffusion 

A  R  Mitchell  and  P  John-Charles,  University  of  Dundee. 

The  non  linear  reaction  diffusion  equation 

IF  '  Ji  (u  +  u(l  "  u) 

and  discretisations  of  it  in  space  and  time  model  many  problems  in 
mathematical  biology.  It  is  shown  by  numerical  experimentation  that 
stable  periodic  travelling  wave  solutions  exist  beyond  the  von  Neumann 
linearised  stability  limit,  even  when  the  grid  sizes  in  space  and  time 
are  small. 


A  Modified  Moving  Finite  Element  Method  for  Moving  Boundary  Problems 
R  0  Moody,  University  of  Reading. 

Element 

A  modified  Moving  Finite/Method  is  described  and  applied  to  three  standard 
one-dimensional  moving  boundary  problems:  one-  and  two-phase  Stefan  and 
the  Crank-Gupta  oxygen  diffusion  with  absorption.  Very  favourable 
results  are  obtained  for  both  the  dependent  variable  and  the  boundary/ 
interface  position,  using  a  small  number  of  elements. 
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Quadrature  Formulae  for  Cauchy  Principal  Value  Integrals  of  Oscillatory 
Kind 

G  E  Okecha,  University  of  Calabar,  Nigeria. 

The  problem  considered  is  that  of  evaluating  nisnerically  an  integral  of 
fl  iwx 

the  form  e  f(x)  dx,  where  f  has  one  simple  pole  in  the  interval 
-1 

1-1,1].  Modified  forms  of  the  Lagrangian  interpolation  formula,  taking 
account  of  the  simple  pole  are  obtained,  and  form  the  bases  for  the 
numerical  quadrature  rules  obtained.  Further  modification  to  deal  with 
the  case  when  an  abscissa  in  the  interpolation  formula  is  coincident  with 
the  pole  is  also  considered.  An  error  bound  is  provided  and  some  numerical 
examples  are  given  to  illustrate  the  formulae  developed. 


On  Certain  Minimal  Polynomials  Related  to  the  Iterative  Solution  of 
Linear  Equations 

G  Opfer,  University  of  Hamburg,  F.R.  Germany. 

The  iteration  scheme  x:+p  «  Xi  -  a: (Ax :  -  b) ,  i  •  1,2,...,  for  solving 
the  (real  or  complex)  linear  N  x  N  system  Ax  -  b  is  particularly  well 
suited  for  parallel  computation,  since  all  components  of  *j+p  can  be 
computed  independently  and  therefore  simultaneously. 

He  assume  that  a  set  S  which  contain  the  spectrum  o(A)  of  A  is 
known.  The  above  iteration  converges  if  the  parameters  dj  are  chosen 

such  that  II  P II  oo  <  wher«  p(*)  "  Rj„j.  d  “  Oj*),  akn+j  ”  a  j »  j"l,2,...,n 

for  any  k  e  W  and  n  «  W  fixed  (cyclic  iteration  with  cycle  length  n) . 

Therefore  the  knowledge  of  the  minimal  polynomials  which  minimise  ||  ’Ho, 
within  all  p  «  Iln  with  p(0)  «  1  is  essential  for  finding  the  above 
mentioned  a's.  In  all  cases  the  maximum  norm  ||  "Ho,  has  to  be  computed 
on  S. 

For  n  «  l  and  arbitrary  sets  S  the  problem  has  been  solved  by 
F.  SCHOBF.R  and  the  lecturer  (LAA  58(1984)  ,343-361) .  We  shall  report  on 
results  mainly  for  n*2  and  rectangular  inclusion  sets  S.  These  sets 
are  particularly  well  suited  for  eigenvalue  inclusions  due  to  a  theorem 
by  BENOIXSON. 
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A  Mew  AjVy^ach  to  Saaah's  Parallel  Eigenvalue  Algorit.—' 
M  H  C  Paardekooper,  Tilburg  University,  The  Netherlands. 


This  paper  presents  an  improved  version  of  Saneh's  parallel  eigenvalue 
algorithm.  It  constructs  a  sequence  {A^}  with  Ak*l“®it  *kSk>  ken 
where  the  normreducing  transformation  matrix  Sw  is  a  direct  sum  of 
uniisodular  two-dimensional  shears.  These  Jacobi-like  shear  transforms 
can  be  carried  out  concurrently.  An  appropriate  choice  of  the 
non-trivial  elements  of  together  with  a  well-chosen  pivot-strategy 

provides  that  {A^}  tends  to  normality:  A^aJ-A^A^-*-  0  (  k-*-  °° ) .  As  a 

consequence  of  the  invariance  of  the  Frobenius  matrix  norm  under  orthogonal 
transformation  the  norm-reduction  should  be  discussed  in  theoretical  terms 
that  are  invariant  under  these  transformations.  In  our  modification  of 
Sameh's  algorithm  the  problem  formulation  with  so-called  Euclidean  parameters 
improves  and  simplifies  the  procedure  and  the  formulae  in  the  related 
minimization  problem. 

Different  aspects  of  the  method  (global  convergence,  speed  of  convergence) 
and  its  implementation  (communication,  storage  strategy)  are  discussed. 


Gaussian  Quadrature  Formulas  for  Splines  with  an  Application  to 
Integral  Equations 

L  Pretonus*;  D  F  Laurie;  D  Eyre,  NRIMS,  South  Africa. 

The  existence  of  Gaussian  quadrature  formulas  for  splines  was  shown,  for 
example,  by  [Karlin  and  Studden,  1966)  and  also  [Hicchelli  and  Pinkus, 
1977).  However,  the  computational  construction  of  these  formulas  received 
little  attention  in  the  literature.  This  may  be  a  consequence  of  the 
fact  that  aust  methods  for  constructing  the  classical  Gaussian  quadrature 
formulas  rely  heavily  on  the  powerful  theory  of  orthogonal  polynomials 
and  therefore  do  not  apply  in  the  case  of  splines.  In  this  presentation 
we  discuss  the  construction  of  Gaussian  quadrature  formulas  for  splines 
and  then  successfully  apply  them  to  integral  equations. 

S.  Karlin  and  W.  J.  Studden,  1966.  Tchebycheff  systems:  with 
applications  in  analysis  and  statistics.  New  York:  Interscience 
Publishers,  a  division  of  John  Wiley. 

C.  A. Hicchelli  and  A.  Pinkus,  1977.  Moment  theory  for  weak  Chebychev 
systems  with  applications  to  monosplines,  quadrature  formulae  and  best 
one-sided  L^-approxiimtions  by  spline  functions  with  fixed  knots. 

SIAM  Journal  of  Mathematical  Analysis,  8,  206-230. 
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Solving  Nonlinc _ 'Equality  Constrained  OptimisationProbleafr  iing  the 

Truncated  Newton  Method  with  Automatic  Differentiation 
R  C  1’rice,  Hatfield  Polytechnic. 

In  this  paper  we  are  concerned  with  nonlinear  equality  constrained 
optimisation  problems.  The  minimum  of  the  Di  Pillo-Grippo  penalty 
function  occurs  at  the  minimum  of  the  corresponding  constrained  problem. 
In  previous  work  we  have  applied  the  truncated  Newton  Method  to 
unconstrained  problems  [1],  where  the  required  derivatives  have  been 
obtained  by  automatic  differentiation  [2].  In  this  paper  we  extend 
these  ideas  and  apply  them  to  the  Di  Pillo-Grippo  penalty  function. 

All  the  derivatives  are  obtained  in  such  a  way  that  the  hessians  of  each 
constraint  are  stored  one  at  a  time  and  not  simultaneously,  thus  making 
the  method  low  on  storage  requirements. 

Encouraging  results  tiave  been  obtained  when  the  method  has  been  applied 
to  a  series  of  test  problems. 

[1]  Numerical  experience  with  the  truncated  Newton  method.  Technical 
Report  No.  169,  NOC,  The  Hatfield  Polytechnic.  (To  appear  in 
Journal  of  Optimisation  Theory  and  Applications). 

[2J  The  truncated  Newton  method  for  sparse  unconstrained  optimisation 
using  automatic  differentiation.  Technical  Report  No. 170,  NOC, 

The  Hatfield  Polytechnic.  (To  appear  in  Journal  of  Optimisation 
Theory  and  Applications). 


An  Alternative  to  Chan's  Deflation  for  Bordered  Systems 
J  D  Pryce,  University  of  Bristol. 


Linear  systems  like 


vl  .  l"A  bl  N 
M  LyJ  '  Lc  dJ  LyJ 


frequently  occur  in  numerical  continuation  and  other  contexts.  Here 
A  is  nxn,  G  is  n  *1,  c  is  l  *n  and  d  is  a  scalar,  and  A  is 
structured  (sparse,  banded,  etc.)  while  the  'border'  b,  c  and  d  is 
usually  full.  The  complete  (n+l)  »  (n+1)  matrix  M  is  assumed  to  be 
well  conditioned  but  A  may  be  arbitrarily  near  to  singularity. 

However,  for  efficiency's  sake  we  wish  to  use  a  solver  for  linear  systems 
of  the  structure  of  A,  as  a  black  box.  How  to  do  this  and  come  out 
with  accurate  answers  even  when  A  is  ill-conditioned? 

The  talk  will  present  an  alternative  to  the  frequently  used  deflation 
method  of  T.  F.  Chau.  The  new  method,  due  to  W.  Govaerts  and  the 
speaker,  is  simpler  to  program  and  in  a  large  number  of  tests  on  moderate¬ 
sized  matrices  has  not  been  found  to  fail. 


Work  on  the  case  of  borders  of  width  k,  1  <  k  <<  n,  will  also  be  mentioned. 
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A  Numeri1,  ,  Stratified  Model  for  Hurricane  Generated  rrents 
K  K  Puri*;  C  Cooper,  University  of  Maine,  U.S.A. 

A  multi-layer  model  is  considered  to  study  hurricane  generated  currents 
in  an  ocean  environment.  It  consists  of  four  layers.  The  upper  most 
mixed  layer  admits  the  presence  of  a  free  surface  and  hence  of  the 
atmospheric  pressure  gradients.  Its  interaction  with  the  lower  layer 
permits  its  thickening  as  a  result  of  entrainment  due  to  the  shear 
instabilities. 

The  time  dependent  equations,  appropriate  to  the  four  layers  which  take 
into  account  the  Coriolis  effect  are  solved  numerically  using  an  explicit 
finite  difference  scheme.  The  results  are  compared  with  existing 
current  models  which  this  model  subsumes.  Also,  barotropic  and  baroclinic 
responses  are  sought,  to  the  presence  of  hurricanes  Fredricks  and  Allen. 


A  Spline  Approximation  He t hod  for  Solving  System  of  O.D.E. 
I  Refat,  Ministry  of  Education,  Kuwait.  ’ 


In  this  paper  we  consider  the  system  of  second  order  ordinary  differential 
equations  with  initial  conditions.  We  constructed  a  spline  function, 
which  is  not  necessarily  polynomial  spline,  approximating  the  solution  of 
this  problem  and  giving  optimal  error  bounds. 


Accuracy  of  Difference  Schemes  for  Solving  the  Wave  Equation 
R  Renaue,  Arizona  State  University,  U.S.A. 


Time  stepping  schemes  for  the  solution  of  the  wave  equation  are  analysed. 
It  is  demonstrated  that  the  determination  of  stability  and  order  of 
accuracy  properties  of  these  schemes  can  be  posed  as  a  rational 
approximation  problem.  Applying  order  star  techniques  a  bound  on 
accuracy  is  derived.  Further,  all  of  the  schemes  attaining  this  bound 
are  stable.  Similar  analysis  leads  to  a  bound  on  accuracy  of  full 
discretisations  of  the  wave  equation.  The  explicit  methods  are  stable 
for  Courant  numbers  less  than  one. 


Choices  in  Contractivi ty  Theory  of  One-Leg  and  Multistep  Methods 
J  Sand  ,  University  of  Copenhagen,  Denmark. 

Contractivity  theory  has  in  the  past  been  used  as  a  rather  powerful 
tool  for  establishing  stability  results  for  the  solution  of  ordinary 
differential  equations  (ODE's).  The  applicability  of  the  tool  is, 
however,  strongly  dependent  on  three  choices,  which  have  to  be  made, 
viz. 

-  The  class  of  ODE's  to  be  considered. 

-  The  class  of  methods  to  be  applied. 

-  The  type  of  norm  in  which  contractivity  is  to  hold. 

Restricting  ourselves  to  the  classes  of  one-leg  and  linear  multistep 
methods,  we  shall  review  some  of  the  choices  which  have  been  examined 
in  the  past  and  suggest  alternatives  which  in  certain  cases  lead  to 
less  pessimistic  results. 


) 

Stability  Phenomena  Arising  in  the  Numerical  Study  of  a  Fluidized  Bed 
J  M  Sanz-Serna*;  L  Abia,  Universidad  de  Valladolid,  Spain. 

When  numerically  studying  an  evolutionary  partial  differential  equation 
modelling  a  fluidized  bed,  I.  Christie  and  6.  Ganser  found  a  number  of 
remarkable  stability  phenomena  (see  their  contribution  to  this  conference). 
Crank-Nicolson  /implicit  midpoint  time  discretizations  of  reasonable 
finite  difference  or  Galerkin  space  discretizations  turned  out  to  be 
unstable,  while  schemes  based  on  the  backward  difference  formulae  behaved 
in  a  stable  manner  only  if  the  time  step  was  chosen  to  be  lafge  relative 
to  the  mesh  size  in  space. 

In  the  present  contribution  we  analyze  in  detail  these  phenomena  and 
suggest  new  schemes  which  do  not  suffer  from  stability  problems.  The 
nonlinear  stability  and  convergence  of  the  new  methods  is  studied  within 
a  framework  developed  by  J.  C.  Lopez-Harcos  and  one  of  the  present  authors 
(J.M.S.). 


Best  Parametric  Interpolation  by  Cubic  Splines  in  IR^ 

K  Scherer,  University  of  Bonn,  W.  Germany. 

The  problem  of  finding  a  curve  by  passing  through  prescribed  points  in 
lRd  i3  considered  such  that  the  kinetic  energy  of  a  particle  moving 
along  this  curve  is  minimized.  Since  the  moment  of  passing  is  free 
(within  a  given  time)  this  results  in  a  nonlinear  problem.  The 
uniqueness  of  its  solution  is  proved  by  showing  that  it  is  equivalent 
to  a  convex  programming  problem.  Also  the  shape  preserving  properties 
I  of  the  solution  are  discussed. 


High  Wave  Modulations  in  Numerical  Solutions  of  the  Korteweg-De  Vries 
and  Modified  Korteweg-De  Vries  Equations 

S  U  Schoombie,  University  of  the  Orange  Free  State,  South  Africa. 

The  results  of  numerical  experiments  are  presented,  showing  a  solitary 
sawtoothed  wave  packet  emerging  from  rectangular  pulse  initial  data. 

In  the  case  of  the  MKdV  -  equation  an  analysis  shows  such  wave  packets 
to  have  envelopes  satisfying  a  different  MKdV  -  equation.  Thus 
sawtoothed  envelope  solitons  become  a  possibility.  These  phenomena 
are  due  to  the  nature  of  the  space  discretization  used,  which  is  a 
central  difference  scheme.  A  multiple  scales  analysis  shows  that 
modulated  sawtoothed  waves  of  small  amplitude  satisfy  a  linear  dispersive 
difference  equation. 
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EMIT  -  Ai  elicit  Runge-Kutta  Hyatrom  Integrator  tor  md  Order 
Initial  Vf.ue  Problems 

F  W  Sharp*;  J  M  Fine,  University  of  Toronto,  Canada. 

In  this  paper  we  discuss  an  integrator  we  have  developed  for  second 
order  initial  value  problems  for  the  fora  y"  ”  f(x,y),  y(xg)  “  y„, 
y* (xQ)  -  y'g.  The  integration  is  performed  using  an  explicit  (4,5) 
Runge-Kutta Nystrom  pair  with  error  control  on  y  and  yf.  Local 
extrapolation  is  enq>loyed  and  y  and  y'  are  estimated  at  the  output 
points  using  interpolants  that  do  not  require  extra  derivative 
evaluations.  We  begin  by  summarizing  properties  of  the  (4,5)  pair 
and  the  interpolants.  Following  this  the  performance  of  the  integrator 
is  illustrated  for  several  test  problems  on  a  scalar  machine.  Finally 

we  discuss  how  the  integrator  can  be  modified  to  take  advantage  of  the 
vectorizing  and  sxiltitasking  capabilities  of  a  CRAY  X-MP/2. 


Solving  Linear  Partial  Differential  Equations  by  Exponential  Splitting  (I) 

Qin  Sheng,  University  of  Cambridge. 

Let  A,,  A.,  ....  Ay  be  square  matrices  which  do  not  coumute.  We 
consider  approximations  to  the  matrix  exponential  H(t)  -  exp ( t (A  +A  +. . .A  )  } 
k  1  2  N 

of  the  form  E  y^E^,  w*iere  each  y.  i»  ■  multiplying  factor,  and  each 
k“l  K 

\  i*  a  product  of  terms,  each  term  having  the  form  exp(atA  }  for  some 
a  >  0  and  1  S  n  £  H.  This  form  is  relevant  to  semi-discreSization 
methods  for  linear  partial  differential  equation  problems,  and  it  produces 
systems  which  are  easy  to  solve.  The  accuracy  and  stability  of  this 
splitting  approximation  are  studied.  It  is  shown  that,  even  if  the 
number  of  terms  and  the  number  of  values  of  k  are  chosen  to  be  large, 

Che  highest  order  of  a  stable  approximation  is  two.  Numerical  examples 
are  given. 


Upwinding  on  Triangular  Elements 
S  Sigurdsson,  University  of  Iceland. 

Upwinding  techniques  for  one-dimensional,  steady-state,  convective 
transport  equations  in  the  framework  of  Petrov-Calerkin  methods  have 
been  studied  extensively  by  many  authors.  However,  the  question  of 
extending  these  techniques  to  two-dimensional  problems,  in  particular 
in  the  case  of  irregular  triangular  finite  elements  seems  as  yet  not 
to  have  been  fully  resolved.  We  present  different  viewpoints  of  this 
question  with  special  enq>hasis  on  linear  basis  functions  and  describe 
techniques  that  have  been  used  for  solving  groundwater  contamination 
problems. 


»• 
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A  l.ocally  Mass  Conserving.  Quadratic  Velocity.  Linear  Pressure  Element 
for  Incompressible  Flow 
IJ  J  Silvester,  IIMIST. 

One  of  the  most  popular  triangular  elements  for  2D  incompressible  flow 
problems  is  the  continuous  quadratic  velocity,  continuous  linear 
pressure  element.  The  principal  attraction  of  the  element  is  its 
uniform  stability  coupled  with  its  cost-effectiveness.  Despite  the 
element’s  general  acceptability,  serious  doubts  concerning  its  perforaiance 
have  also  been  expressed.  In  particular,  the  inherent  lack  of  continuity 
concomitant  with  a  continuous  pressure  approximation  is  known  to  give 
rise  to  non-physical  solutions  if  the  flow  problem  is  sufficiently 
'difficult'.  For  example,  poor  results  using  the  element  to  model 
recirculating  isothermal  flows  and  strongly  coupled  thermal  flows  have 
been  presented  previously.  A  less  well-known  class  of  problems  which 
t  also  demonstrate  this  behaviour  is  that  of  non-Newtonian  flow. 

To  improve  the  approximation  of  such  flows  a  'hybrid'  element  having  two 
superposed  pressure  fields  has  been  analysed  and  tested  numerically. 

In  particular,  supplementing  the  basic  pressure  field  with  a  piecewise 
constant  pressure  field  has  been  shown  to  lead  to  an  improved  mixed 
approximation  which  is  still  uniformly  stable.  By  way  of  illustration, 
the  performance  of  the  modified  element  is  contrasted  with  that  of  the 
conventional  element  for  a  rudimentary  non-Newtonian  swirling  flow. 


Piecewise  Polynomials  and  Integrodif ferential  Equations 
1  M  Snyman*; S  A  Sofianos,  University  of  South  Africa. 

The  integrodif ferential  equation  derived  in  the  framework  of  the 
llyperspherical  Harmonics  method  to  describe  the  bound  state  of  a 
three-boson  system  in  configuration  space,  has  been  solved  using 
B-splines.  The  equation  includes  the  strong  effects  of  the  hypercentral 
part  of  the  two-body  potential  and  has  apparent  singularities  as  well  as 
6-function  type  behaviour  in  the  asymptotic  region.  B-splines  proved 
to  be  particularly  suitable  to  handle  these  difficulties.  The  Galerkin 
method  is  used  together  with  soft  and  hard  core  two-body  interactions  to 
obtain  upper  and  lower  bounds  for  the  binding  energy  of  the  system. 
Various  numerical  aspects  of  the  problem  such  as  stability  of  the  method, 
choice  of  mesh  points  and  the  order  of  the  B-splines  are  discussed. 


A  Generalisation  of  Codunov'a  Method  via  Randomly  Samp  leu  Intercell  Fluxes 
E  F  Toro,  Cranfield  Institute  of  Technology. 


A  large  variety  of  Modern  computational  techniques  for  quasi-linear 
hyperbolic  equations  make  use  of  the  solution  of  the  Riemann  problem. 

Here,  yet  another  method  is  presented  in  which  to  update  the  solution  at 
cell  i  a  leap-frog  finite  difference  scheme  is  used,  i.e. 


where  intercell  fluxes 
problem  RP(i,  i+1) 


f"+1 

ri*l 


are  evaluated  at  the  solution  of  the  Riemann 


at  a  randomly  given  point  P.  -  (x,  tn+i) 
.  Here  6n  is  a  quasi-random  number  lA 


l)  with 

the  interval 

C-|T,  JT]  and  T  is  an  interval-length  sampling  function  varying  between 
T  ■  0  and  T  ■  1.  Every  value  of  T  in  [0,1]  gives  a  difference  scheme 
in  particular  T  «  0  reproduces  Godunov's  first  order  method,  T  *  1 
gives  a  second  order  method  that  can  be  interpreted  as  a  random  Lax-Uendroff 
Method. 


Suitable  construction  of  the  function  T  leads  to  a  method  that  ranges 
through  the  family  of  schemes  described  above  giving  higher  accuracy  than 
Godunov's  Method  and  which  is  also  monotone. 

Applications  to  the  Euler  equations  in  one  and  two  dimensions  are  made. 


On  the  Convergence  of  the  Product  Approximation  in  Semilinear  Problems 
t  Tourigny,  University  of  Dundee. 


By  aieans  of  a  concrete  example,  we  will  first  explain  what  the  so-called 
product  approximation  in  the  finite-element  method  consists  of.  We 
will  then  review  the  existing  convergence  theory  for  this  technique. 

At  the  moment,  the  theory  covers  such  problems  as:  nonlinear  Klein-Gordon 
equations,  nonlinear  Schrodinger  equations,  nonlinear  parabolic  equations 
and  semilinear  tun-point  boundary  value  problems. 


Is  Gaussian  Elimination  Stable? 

L  N  Trefethen*;  R  S  Schreiber,  MIT,  U.S.A. 

Wilkinson  showed  that  Gaussian  elimination  with  partial  pivoting  is  stable 
so  long  as  the  "growth  factor"  is  not  too  large.  But  he  also  showed  that 
in  the  worst  case  it  is  too  large  -  2n-^  for  a  matrix  of  dimension  n. 

Thus  Gaussian  elimination  is  unstable  in  the  worst  case. 

Why  then  does  this  tried-and-true  algorithm  work  so  well  in  practice? 

Some  people  have  proposed  that  the  matrices  that  arise  in  practical 
problems  must  be  exceptionally  well-behaved.  In  this  talk  we  show  that 
there  is  no  need  for  that  hypothesis.  In  fact,  Gaussian  elimination  is 
stable  on  average,  and  the  average-case  growth  factor  is  no  worse  than 
n.  We  support  these  claims  with  empirical  evidence  and  the  beginnings 
of  a  theoretical  explanation. 


to  Convection-Diffusion  Equations 
J  G  Verwer,  Centre  for  Mathematics  and  Computer  Science,  The  Netherlands. 

The  hopscotch  method  is  a  time  stepping  scheme  applicable  to  wide  classes 
of  spatially  discretized,  multispace  dimensional,  time-dependent  partial 
differential  equations.  In  this  contribution  attention  is  focussed  on 
the  simple  odd-even  hopscotch  method.  Our  aim  is  to  present  some 
interesting,  linear  stability  results  of  this  method  for  convection- 
diffusion  equations,  where  the  space  discretization  is  carried  out  by 
standard  synmetrical  and/or  one-sided  finite  differences.  We  present 
two  sorts  of  results.  The  first  result  is  concerned  with  von  Neumann 
stability  in  the  multi-dimensional  case.  Here  we  derive  expressions  for 
the  critical  time  step  and  show  that  in  certain  cases  an  increase  of 
diffusion  may  render  the  process  unstable,  an  observation  which  is  in 
clear  contrast  to  the  conmon  practice.  This  strange  phenomenon  can 
occur  only  in  the  higher  dimensional  case.  The  second  result  is 
concerned  with  the  spectral  condition,  but  only  in  one  space  dimension. 

We  derive  the  critical  time  step  values  and  conclude,  in  accordance  with 
known  results  for  the  explicit  Euler  rule,  that  the  spectral  condition 
is  misleading  in  the  sense  that  it  does  not  prevent  large  error  growth. 


A  Type-Insensitive  ODE  Solver 

K  C  Wade  ;  C  W  Richards  ;  M~G  Everett  ,  Thames  Polytechnic,  London. 

A  variety  of  efficient  integration  algorithms,  in  the  form  of  subroutine 
libraries,  are  available  for  the  solution  of  initial  value  problems  in 
UUE's.  However  it  is  left  to  users  to  select  the  most  appropriate 
algorithm  for  their  problem.  Frequently  inexperienced  users  select 
unsuitable  algorithms  such  as  attempting  to  integrate  stiff  systems 
using  explicit  Runge-Kutta  methods.  The  situation  is  further  complicated 
when,  as  often  happens  in  engineering  problems,  the  characteristics  of 
the  equations  change  during  the  integration  range. 

This  paper  describes  a  numerical  code  based  on  Runge-Kutta  methods  that 
assesses  the  characteristics  of  the  equation  throughout  the  integration 
range  and  automatically  selects  the  most  appropriate  integration  method. 

A  total  of  three  explicit  Runge-Kutta  and  four  implicit  Runge-Kutta  methods 
are  available  giving  a  variable  step,  variable  order  algorithm.  The 
strategy  employed  in  selecting  order  and  in  switching  beteen  explicit  and 
implicit  methods  is  described.  To  ensure  that  the  switching  strategies 
are  efficient  an  embedding  technique  is  eng>loyed,  which  requires  no 
additional  function  evaluations. 


Numerical  results  are  presented  which  demonstrate  the  effectiveness  of 
the  switching  technique  and  overall  efficiency  of  the  code. 


3ft 


Global  Preconditioning  for  Galerkin  Finite  Element  Equanuns  Using 
Element-by-Element  Analysis 
A  J  Wathen,  University  of  Bristol . 


It  is  traditional  in  Finite  Element  calculations  to  assemble  information 
from  the  individual  elements  into  large  global  matrix  systems,  then  to 
solve  using  elimination  methods.  The  matrix  is  usually  large  and  is 
generally  assumed  to  be  banded  about  the  diagonal  or  to  have  some  sparsity 
structure.  Such  a  general  approach  ignores  the  specific  structure  which 
is  due  to  the  construction  (assemblage)  of  the  matrix  itself.  In  this 
paper  we  show  how  properties  of  the  matrix  which  are  useful  in  the  solution 
process  can  be  derived  simply  from  observations  regarding  the  assemblage 
of  the  global  matrix  from  individual  element  matrices.  Specifically  we 
show  how  to  derive  eigenvalue  bounds  for  symmetric  finite  element 
coefficient  matrices  with  simple  preconditioners  which  are  formed  also  by 
an  assembly  process.  These  bounds  on  the  eigenvalues  of  the  preconditioned 
global  matrix  are  derived  element-by-element  and  have  been  used  by  the 
author  to  guarantee  rapid  convergence  of  the  Preconditioned  Conjugate 
Gradient  Method  for  the  solution  of  finite  element  equations.  We  give 
results  for  diagonal  and  tridiagonal  preconditioning  of  Galerkin  mass 
and  stiffness  matrices  for  many  common  types  of  element  in  1-  2-  and 
3-dimensions  as  well  as  an  example  of  the  wider  use  of  the  technique. 

It  will  also  be  shown  how  this  work  explains  some  of  the  success  of 
Hughes'  element-by-element  algorithm  for  elliptic  problems  (Hughes,  T.J.R., 
Levit,  I.  and  Winget,  J.,  Comput.  Meths.  Appl.  Mech.  Engng.  36,  241-254, 
1983). 


The  Smallest  Perturbation  of  a  Submatrix  which  Lowers  the  Rank  of  a  Matrix 
G  A  Watson,  University  of  Dundee. 


The  problem  is  considered  of  determining  the  smallest  perturbation  of 
a  submatrix  of  a  partitioned  rectangular  matrix  which  achieves  for  that 
matrix  a  given  lower  rank.  For  different  classes  of  matrix  norms,  it 
is  shown  how  the  solution  of  the  problem  can  he  obtained  through  an 
equivalent,  often  simpler  problem. 


Hie  Eigenvalues  of  Spectral  Differentiation  Matrices 
J  A  C  Weideman*;  L  N  Trefethen,  MIT,  U.S.A. 

Spectral  methods  for  partial  differential  equations  are  based  upon  the 
approximation  of  space  derivatives  by  global  trigonometric  or  polynomial 
interpolation  followed  by  differentiation  of  the  interpolant.  This 
process  is  equivalent  to  the  multiplication  of  the  data  vector  by  a 
matrix  D.  The  eigenvalues  of  D  are  of  fundamental  importance  to  the 
stability  and  other  numerical  properties  of  such  methods.  This  talk 
surveys  old  and  new  results  on  the  distribution  of  these  eigenvalues  for 
fundamental  spectral  problems  involving  first-  or  second-order 
differentiation  in  equispaced  or  Chebyshev  points. 
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An  Algorithm  for  Minimizing  a  Quadratic  Function  with  Quadratic  Constraints 
Y  Yuan,  University  of  Cambridge. 


We  study  the  problem  of  minimizing  a  general  quadratic  function  with 
special  quadratic  constraints.  This  problem  is  a  sub-problem  in  some 
trust  region  algorithms  for  general  nonlinear  optimization.  An 
algorithm  is  presented. 


The  Rate  of  Convergence  of  Sliding  Block  Jacobi  Schemes 
Yong  l.i ;  W  Ferguson*,  Southern  Methodist  University,  U.S.A. 

The  sliding  block  Jacobi  schemes  considered  in  this  paper  are  q-line  block 
Jacobi  schemes  in  which  the  last  i  lines  of  one  block  are  overlapped  with 
the  first  r— lines  of  the  next  block.  Analytic  formulas  for  the  spectral 
radius  of  these  sliding  block  Jacobi  schemes  are  presented  for  a  class 
of  model  problems  which  includes  the  classical  Poisson  equation  on  a 
rectangle.  Well  known  formulas  relate  the  spectral  radius  of  the  Jacobi 
scheme  to  the  rate  of  convergence  of  associated  SOR  and  CG  schemes. 
Sliding  block  Jacobi  schemes  may  form  the  basis  of  robust  and  efficient 
iterative  schemes  for  solving  elliptic  partial  differential  equations  on 
parallel  computers  because  their  rate  of  convergence  increases  rapidly 
with  increasing  values  of  q  and  r. 


Iterative  Methods  for  Nonlinear  Operator  Equations  of  Monotone  Type 
You  Zhaoyong*; Xu  Zongben,  Xi'an Jiaotong  University,  P.R.  China. 


Many  problems  arising  naturally  in  the  field  of  differential  equations, 
integral  equations,  variational  methods  and  optimizations  involve  the 
solution  of  equations  associated  with  nonlinear  operators  of  monotone 
type.  In  this  paper,  we  systematically  develop  iterative  methods  for 
solving  such  different  type  of  monotone  operator  equations,  which  includes: 

-  Relaxation  iterative  processes  for  strongly  monotone  operators 

-  A  second  order  difference  regularization  approach  to  maximal  monotone 
operators 

-  A  class  of  iterative  method  with  finite  termination  property 

-  Ergodic  convergence  of  Ishikawa's  method  for  general  monotone  operators. 
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Upper  and  Lower  Bounds  to  Eigenvalues  of  Eigenvalue  Problems  with 
Partial  Differential  Equations 

S  Zimmermann,  Technische  Universitat  Clausthal,  Germany. 


This  paper  deals  with  eigenvalue  problems  M$  *  XN$  (X  e  IR  eigenvalue, 

$  €  H  eigenvector),  where  M,  N  are  linear  syntnetric  operators  in  a  llilbert 
space  H.  Accurate  upper  bounds  to  eigenvalues  can  often  be  obtained  by 
applying  the  Ray leigh-Ritz  method.  For  calculating  corresponding  lower 
bounds  one  can  sometimes  employ  the  Lehmann-Maehly  method,  which  contains 
the  well-known  Temple  formula  as  a  special  case.  A  generalization  of 
this  method  has  been  proposed  by  Goerisch  recently,  it  possesses  a  much 
larger  range  of  application.  The  convergence  of  both  Lehmann  and 
Coerisch's  bounds  is  studied.  Numerical  results  are  given  for  the  plate 
prob lems 


A1#  “  -XA$  in  11, 


0  on  3(1 


(buckling) , 


A2*  -  X*  in  (1, 


0  on  30 


(vibration) . 


